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As a function of applied field, we find a rich variety of ordered and partially-ordered vortex lattice 
configurations in systems with square or triangular arrays of pinning sites. We present formulas 
that predict the matching fields at which commensurate vortex configurations occur and the vortex 
lattice orientation with respect to the pinning lattice. Our results are in excellent agreement with 
recent imaging experiments on square pinning arrays [K. Harada et al., Science 274, 1167 (1996)]. 
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I. INTRODUCTION 

The concept of an elastic lattice interacting with a 
rigid substrate lattice producing commensurate and in- 
commensurate transitions when the periodicities of the 
two lattices match or mismatch is found in numerous 
condensed matter systems, including atoms adsorbed on 
surfaces Q, layered superconductors Q, superconduct- 
ing networks and Josephson-junction arrays ||, colloids 
HI), and magnetic bubble arrays interacting with pat- 
terned substrates ||. Recently, increased interest has 
been focused on superconducting systems with well de- 
fined square or triangular periodic pinning arrays in 
which vortices can be trapped both at individual pinning 
sites and also at the interstitial regions between pinning 
sites [fz| — |l3|| - Due to the interstitial pinning, the vortex 
system differs significantly from systems such as atoms 
on surfaces or Josephson-junction arrays. For instance, 
in the latter systems, the potential substrate has an egg- 
carton form, so the atom or vortex lattice has the same 
configuration whenever it matches the underlying poten- 
tial. Recent direct imaging experiments (?|] and simula- 
tions [ fl2| indicate that ordered vortex configurations in 
samples with periodic pinning can vary at each matching 
field, producing a remarkable variety of stabilized vortex 
lattices which are quite distinct from those found in su- 
perconducting networks Q and other systems. 

Since highly ordered commensurate lattices can be 
more strongly pinned than incommensurate lattices 
|7|,p[~p^[ , a determination of how different matching con- 
figurations affect the overall pinning of the vortex lat- 
tice could be useful for technological applications of su- 
perconductors. For instance, enhanced pinning at cer- 
tain matching fields has been verified with the observa- 
tion of peaks in easily measurable magnetization curves 
§,001, including high-T c materials Q- 

Imaging experiments have so far only probed up to 
the fourth matching field and have only examined square 
pinning arrays . A general characterization of the vor- 
tex matching patterns as a function of arbitrary matching 



densities for square and triangular pinning arrays has not 
been done up to this point. 

We have performed a series of large scale simulated 
annealing as well as flux-gradient-driven [ |l4| molecular 
dynamics simulations of vortices interacting with square 
and triangular arrays of small pinning sites for very high 
fields (up to the 28th matching field), and for a wide 
range of pinning parameters and system sizes. 

Our results show that the vortex lattice (VL) is highly 
ordered only at certain matching fields (MFs) and can 
have various orientations with respect to the underlying 
pinning array. At some MFs the VL is actually disor- 
dered. The enhancements of M(H) are most noticeable 
for fields less than the second matching field; however, we 
find some evidence of small enhancements of M(H) for 
higher fields. Square and triangular arrays produce dif- 
ferent sequences of ordered matching fields at which the 
pinning is enhanced. At some MFs, we find novel vortex 
arrangements with translational order only along certain 
directions. Our numerical results are in excellent agree- 
ment with recent low-field experiments on square pinning 
arrays J?]]. Moreover, using geometrical arguments that 
take into account the constraints of the pinning array, 
we derive simple formulas for the ordered MFs and for 
the orientation of the VL with respect to the square or 
triangular pinning array. 



II. SIMULATION 

We perform simulated annealing simulations for a 2D 
transverse slice (in the x-y plane) of a superconducting 
3D slab by numerically integrating the overdamped equa- 
tion of motion of rigid 2D vortices: 



rpv., 



= f,: = fT 



(1) 



The term fj is the total force per unit length acting on 
vortex i. The force due to the interactions with other 
vortices is 
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*T= E fa* (^hr^) f« (2) 

3=1 ^ ' 

where TV,; is the number of vortices, fy = (r$ — ij)/|rj — 
Tjl, and we take 77 = 1. Ki(r/X) is the modified Bessel 
function, A is the penetration depth, and 

f *° 
Jo 8^ 2 A 3 

The Bessel function decays exponentially for r greater 
than A, so for computational efficiency the interaction 
can be safely cut off at 6A. In thin-film superconductors 
the long-range vortex- vortex interaction decays as 1/r 
unlike in 3D bulk superconductors; however, the excellent 
agreement between our results and experiments in thin 
films (jj indicates that our results are valid for both slabs 
and thin films and are general enough to be applicable 
to other systems with repulsive particles on a periodic 
substrate (e.g., colloids). The pinning force is 




(3) 

where O is the Heaviside step function, f p is the max- 
imum pinning force, N p is the number of pinning sites 

and = (rj — rj^ )/|rj — rjj Temperature is added 
as a stochastic term with properties 

</?(*) >=0 (4) 

and 

< f?(t)fj(t) >= 2i 1 k B T8 lJ 8{t~t). (5) 

To find the vortex ground state, we gradually cool a fixed 
number of randomly moving vortices from a high temper- 
ature to T = 0, simulating the field-cooled experiments 
of Ref . M . To examine vortex mobility and features in 
the magnetization curves as a function of applied field 
H , we use flux-gradicnt-drivcn simulations in which only 
the central 2/3 of the sample contains pinning sites. In 
this case, vortex lines are slowly added to the unpinned 
region and force their way into the pinned region (the ac- 
tual sample). Although the vortex system in flux-driven 
simulations is in a non-equilibrium state, almost all of the 
vortex states found by simulated annealing also appear 
in the flux-driven case in parts of the sample. 

We measure all lengths in units of A and fields in 
$o/A 2 , and consider systems from 36A x 36A up to 
72A x 72A in size. The pinning is placed in square or 
triangular arrays at densities between n p = 0.072/A 2 and 
0.81A 2 . The pinning radius is fixed at r p = 0.35A. Pin- 
ning sites this size and smaller trap only one vortex per 
pinning site, which is similar to the experimental situa- 
tion in Ref. M. We consider pinning forces f p varying 
from 0.2/o to fo and examine the VL ordering up to the 
28th MF. 



III. VORTEX LATTICE GROUND STATES: 
SQUARE CASE 

A. Patterns Experimentally Observed 

In Fig. 1 we show a series of VL orderings after anneal- 
ing from our simulations for a square pinning array with 
B = 0.17$ /A 2 for each integer MF up to the 9th MF. 
In (a), at the first MF, all the vortices are trapped at 
the pinning sites so that the overall VL is square. At the 
second MF (b) the interstitial vortices occupy the regions 
in between the pinning sites, so the overall VL is square 
but rotated 45° with respect to the pinning array. At 
the 3rd MF (c), the VL is still highly ordered, with pairs 
of interstitial vortices alternating in position. In (d), at 
the fourth MF, a VL with triangular ordering is observed. 
These structures for the first four MFs correspond exactly 
to those found in direct imaging experiments . We also 
observe ordered VLs at the sub-MFs {B/B$ — 1/4 and 
1/2, where B^ is the vortex density at the first MF), and 
partially ordered VLs at fractional MFs {BjB$ = 3/2 
and B/Bfj, = 5/2) in agreement with experiment (7J. We 
find that the general features of the observed VL config- 
urations up to the fourth matching field are robust for a 
wide range of parameters with 0.2/o < f p < /o, and also 
for 0.072 $ /A 2 < B^ < 0.81 $ /A 2 , for system sizes up 
to 72A x 72A. 

B. Patterns Not Yet Experimentally Observed 

In Figs, l(e-i) we show vortex configurations from our 
simulations that have not yet been observed experimen- 
tally. In (e), at the fifth MF, the overall VL is again 
square and rotated 27° with respect to the pinning lat- 
tice. In (f) a very unusual VL is observed; although 
the VL is neither square nor triangular some ordering is 
still visible. Along the (—1, 1) direction the vortices are 
spaced periodically while in other directions apparently 
periodic distortions can be clearly seen. At the 7th MF, 
in (g), the VL is disordered. In (h), at the 8th MF, the 
VL is nearly triangular. In (i), at the 9th MF, a distorted 
square VL with two different orientations appears, sepa- 
rated by a twin boundary in the middle of the figure. For 
similar systems with lower pin density we have studied up 
to the 28th MF. We see the same VLs already described 
as well as ordered VLs at the 12th and 15th MFs, while 
the vortex configurations at the other MFs have no par- 
ticular ordering. For B/B^ > 15, at high MFs with no 
overall lattice order, the VL contains ordered domains 
separated by grain boundaries of defects similar to those 
observed in Ref. ||]. 
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FIG. 1. Vortex ground states obtained from simulated annealing for a square pinning array with B = 0.17 $o/A 2 , r p = 0.35A, 
and f p /f = 0.625 showing a 12A x 12A subset of a 36A x 36A sample. The flux density is B/B^ = 1 in (a), 2 (b), 3 (c), 4 (d), 
5 (e), 6 (f), 7 (g), 8 (h), and 9 (i). 



Due to numerical constraints we could only look at pin- 
ning densities up to n p = 0.35<f>o/A 2 for 5 < B/B^ < 12, 
and n p = 0.072A 2 for matching fields 12 < B/B^ < 28. 
The vortex patterns observed here are robust for sys- 
tem sizes up to 72A x 72A. The fact that the same pat- 
terns appear for different-sized systems indicates that the 
patterns arise due to commensurability with the pinning 
lattice rather than commensurability with the periodic 
boundary conditions. We should point out that since we 
cannot do infinite-size systems, we cannot conclusively 



rule out finite size effects on the vortex patterns observed. 
Also, in an experimental sample, edges and line/planar 
defects might distort an otherwise periodic VL and cre- 
ate ordered domains that do not extend over the entire 
sample. 
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FIG. 2. Vortex ground states obtained from simulated annealing for a triangular pinning array with B = 0.17 $o/A 2 , 
r p = 0.35A, and f p /f = 0.625/ for a 12A x 12A subset of a 36A x 36A sample. The flux density is B/B^ = 1 in (a), 2 (b), 3, 
(c), 4, (d), 5, (e), 6 (f), 7 (g), 8 (h), and 9 (i). 



IV. VORTEX LATTICE GROUND STATES: 
TRIANGULAR CASE 

In Fig. 2 we show VL configurations up to the 9th 
MF for a system with the same parameters as in Fig. 1 
but with a triangular array of pins. In (a) at the first 
MF, all the vortices are pinned in a triangular lattice. In 
(b) a honeycomb VL forms, with some defects present. 
Ordered triangular VLs are seen at the 3rd (c) and 4th 
(d) MFs with the VL rotated 30° and 0° respectively in 



relation to the pinning array. No clear ordering is seen at 
the 5th (c), 6th (f), and 8th (h) MFs. A triangular VL 
rotated 18° with respect to the pinning array is seen at 
the 7th (g) MF. In (i) at the 9th MF the VL is triangular 
and not rotated with respect to the pinning array. As in 
the samples with square pinning, we observe ordered VLs 
at certain sub-MFs but the VL is disordered at other non- 
MFs. The matching vortex configurations are very robust 
for all the parameters we have investigated, except for 
the honeycomb VL which disappears for weak pinning, 
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fp < 0.3/ . For different samples (n p = 0.072 $ /A 2 , 
f p = 0.625/o) we have studied up to the 28th MF, and 
find ordered triangular VLs at the MFs of order 12, 13, 
16, 19, 21, 25 and 28. The vortex patterns observed for 
the triangular pinning array are robust for a similar set 
of parameters as the square pinning array discussed in 
the previous section. 

V. MATCHING CONDITIONS 

To derive formulas for the fields at which VLs will be 
ordered, for a triangular array of pins, first take a match- 
ing field N and consider any two pinning sites along a 
symmetry axis that have N vortices between them. Now 
find a third pinning site that forms a 60° angle with the 
original two pinning sites, so all three vertices form an 
equilateral triangle. The distances between the third site 
and each of the other two are integer multiples of the 
pinning lattice constant a where 




We label these integers n and to. The distance between 
vortices ojv at a field B is 

a N = 1-075^^. (7) 

By using the law of cosines the distances must obey 

(Na N ) 2 = (ma) 2 + (no) 2 + 2mna 2 cos(60°). (8) 

At the MFs, B = NB,p, so that a N = a/y/N. Substitut- 
ing this into Eq. (8) gives 

N = m 2 +n 2 + nm . (9) 

This equality predicts that for a triangular array of 
pins, an ordered VL will form at values ./V = 1,3,4,7, 
and 9, exactly as seen in our simulations in Figs. 2(a- 
i). Equation (9) also predicts the higher MFs (N — 
12,13,16,19,21,25 and 28) that we have observed nu- 
merically. The honeycomb VL seen in Fig. 2(b) is not 
predicted by Eq. (9) because Eq. (9) only predicts when 
triangular VLs occur. From Eq. (1) we find that the 
angle the VL makes with respect to the pinning array is 

/ V3m/2\ 

6 = arctan '— . (10) 

\n+l/2) V ' 

This equation indicates that for 

B = n 2 B , (11) 

when to = 0, the VL is not rotated with respect to the 
array of pins. This is observed for N = 1, 4, and 9 [shown 



in Figs. 2(a), (d), and (i)], as well as for N — 16 and 
25. For MFs N = 3, where n = 1, to = 1, and N = 
7, where n = 1, to = 2, Eq. (9) gives 9 = 30° and 
19.11° respectively, in agreement with the VLs shown in 
Figs. 2(c,g). We have found that Eq. (10) is valid at least 
up to the 28th MF studied in our simulations. 

We can derive similar conditions for the square pinning 
array, predicting that ordered VLs appear at the Nth MF 
when 

N = m 2 + n 2 (12) 

This equation predicts square VLs for N = 1, 2, 4, 5, 8 
and 9. Indeed, ordered VLs are seen at these fields (see 
Fig. 1). However, only N = 1, 2, and 5 are square in the 
simulation. Moreover, the angle of the VL with respect 
to the array of pins is in principle expected to be 

/ 771 \ 

6 = arctan — (13) 
V n J 

These matching conditions do not always predict the 
right VL ordering observed in simulations. For instance, 
the VLs seen at higher fields N > 9 have triangular 
or distorted square rather than square ordering. These 
equations fail when the VL tendency to remain triangular 
dominates the tendency of the pin array to force a square 
ordering on the VL. This is particularly clear for higher 
fields, N > 9, when the many interstitial vortices are free 
to minimize their energy by forming triangular lattices. 
Equation (9) for the triangular array of pins does not 
have this limitation because both the sample and the VL 
favor a triangular order. 

As we have seen from the simulations, the vortex lat- 
tice is ordered at the matching fields where the commen- 
surability conditions outlined above are met and gener- 
ally disordered where they are not. Several low matching 
fields where these conditions are not met still produce or- 
dered or partially ordered lattices such as the honeycomb 
lattice at the second matching field for the triangular pin- 
ning lattice and the alternating interstitial lattice at the 
3rd field for the square array. This ordering at fields 
not met by our commcnsurability conditions may occur 
due to the pinning being more dominant at lower fields 
so that ordering can be imposed on the interstitial vor- 
tices. For higher fields the vortex configuration for fields 
where commensurability conditions are not met is disor- 
dered or partially disordered. At higher fields B > QB^ 
the vortex-vortex interactions dominate. Here a trian- 
gular vortex lattice is always preferred, so any alternate 
ordering imposed by the pinning does not occur. 
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FIG. 3. Magnetization M(H) and the fraction of vortices with coordination number k, Pk(H), for triangular (a-d), square 
(e-h), and random (i-1), arrays of pins with the same pinning parameters as found in Figs. 1 and 2. In (c) and (g), P5 (P7) 
is represented with a dotted (solid) line. In (a) peaks in M(H) can be seen at the 3rd, 4th, 7th and 9th MFs, which coincide 
with peaks in Pe(H) shown in (b). For the square pinning, small peaks are seen in M(H), (e), at the 2nd, 4th, 5th and 8th 
MFs. The 4th and 5th (small) peaks merge near the 5th MF. The peaks in Pe, (f), indicate that the VL is triangular at the 
4th, 6th and 8th MFs, while peaks in P4, (h), indicate that the VL is square at the 2nd and 5th MFs. Sometimes, the peak is 
slightly shifted (e.g. to 2.5 instead of 2 or 4.5 instead of 4) because of the gradient in the fields. In (g) P5 is the upper curve, 
while in (c) and (k) P5 and P7 follow each other. For the random pinning array at low fields the maximum value \M(H)\ is 
0.0095 <l>o/A 2 , about 1.5 times less than the triangular or square pinning arrays (a,b). For H > 0.4$o/A 2 , the magnetization 
M(H) falls off smoothly while P3 slowly increases as vortex- vortex interactions dominate at higher fields. 
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VI. TOPOLOGICAL ORDER AND 
MAGNETIZATION 

In order to relate the MF configurations to vortex mo- 
bility as well as to experimentally measurable bulk quan- 
tities we present the magnetization M(H) obtained from 
flux-gradient-driven simulations ]T^ | of samples with the 
same pinning parameters used in Fig. 2. In Fig. 3 we plot 
Pk (H) , the fraction of vortices with coordination number 
k obtained from the Voronoi construction, and M(H), a 
useful and common measurement of the net critical cur- 
rent. For a perfect triangular lattice Pq = 1, so any de- 
parture from this indicates a defective lattice. A peak in 
M(H) indicates enhanced pinning. In Fig. 3(a) we find 
that M(H) is very large for fields less than about 2£?^ and 
falls off very rapidly after this. Peaks in both M{H) and 
P 6 {H) appear [Figs. 3(a,b)] at the MFs N = 3,4,7, and 
9 that produced triangular VLs in the simulated anneal- 
ing [Fig. (2)]. From Fig. 2(a), the vortex lattice would 
be expected to form a triangular lattice with the pinning 
substrate at the first matching field, and P 6 would be 
expected to be equal to one. In the flux-gradient-driven 
case shown in Fig. 3(a,b), no peak in M(H) or P§ near 
the first matching field is observed. This is due to the fact 
that the large flux gradient at low fields strongly distorts 
the VL. In the field cooled situation shown in Fig. 2(a) 
there is no gradient in the vortex density to interfere with 
the vortex lattice ordering. For the flux-gradient-driven 
case at higher fields, B > 2B<p, the gradient flattens so 
that the vortex lattice can become commensurate with 
the pinning substrate over a large area. We note that 
even for high matching fields a small flux-gradient will 
always be present so that Pq will always be less than 
one as seen in Fig. 3(b). For weaker pinning, f p < 0.3 fo, 
the flux gradient is reduced at low fields so that a peak in 
M (H) and P 6 can be observed at B/B^ = 1 |§ . We find 
that this behavior is independent of system size. At the 
MFs N = 5 and N = 6 the VL is highly defective with 
Pq{H) dropping as low as 0.5. No peak appears in M (H) 
for the 2nd MF. For systems in which we have studied up 
to the 28th MF we also see some enhancements in M(H) 
and P 6 (H) at the MFs predicted by Eq. (9), although 
the features are washed out at high fields. 

In Fig. 3(e) we show M(H) for a square pinning array 
with the same parameters used in Fig. 1. Again M(H) 
is large for low fields and rapidly falls off after the 2nd 
MF. We can see a dip after the 3rd MF and an over- 
all enhancement in M{H) at the 2nd, 4th, 5th, and 8th 
MFs, although no clear enhancement is seen at the 6th 
and 9th MFs even though the VLs observed through sim- 
ulated annealing at these fields also appear in this flux- 
gradient driven simulation. To examine the evolution of 
the vortices with four nearest-neighbors we consider a 
slightly modified Voronoi algorithm in which the lengths 
of each side of a Voronoi cell are compared. If the length 



of any side is less than one-fourth of the average lengths 
of the other sides, then it is ignored. Pa{H) first shows 
a peak at the second MF when the vortices form the lat- 
tice shown in Fig. 1(b). There is no peak in P±{H) at the 
first MF due to the large flux gradient. Pq(H) shows a 
large peak at the 4th MF that corresponds to the trian- 
gular VL seen in Fig. 1(d). P&{H) then drops rapidly and 
P±{H) increases as the VL gains the square ordering seen 
in Fig. 1(e). Pq(H) rises at the 6th MF and peaks at the 
8th. In square pinning arrays, where we have gone up to 
the 28th MF, small enhancements of M(H) are observed 
for most of the MFs that produced ordered VLs. The 
results indicate that, without directly imaging the VL, it 
could be experimentally possible to deduce the existence 
of the ordered vortex arrays seen here, by looking for a 
specific sequence of peaks in M(H), at least up to the 
5th matching field. Beyond the fifth matching field we 
observe only very small peaks in M(H), which may make 
them difficult to see experimentally. 

Our results are only valid for pins small enough that 
only one vortex can be trapped in each pinning site. 
With triangular pinning, peaks in M(H) should in prin- 
ciple occur for MFs N that satisfy Eq. (9). For square 
pinning arrays, we observe that peaks in M(H) occur 
for MFs given by N = n 2 + m 2 , when N < 10, and 
by N = n 2 + m 2 — 1, when N > 10. This pattern of 
peaks differs from those already seen experimentally us- 
ing periodic pinning arrays with large pinning radii, as 
first shown in Ref. [JL2|. In experiments, peaks in M(H) 
are usually observed at every MF due to multiple vortices 
being trapped in pinning sites |^|,[ll) ■ 

To compare the effects of random pinning to square 
and triangular arrays, in Fig. 3(i) we plot M(H), and 
in Fig. 3(j-l) we plot Pk(H) for a sample with the same 
pinning parameters as in Fig. 3(a,e) except the pinning 
sites are placed in a random array. It can clearly be seen 
that most of the peaks in Pk and M(H) are washed away 
with M(H) having a smooth falloff after the peak and the 
strong variations in Pk lost. The fraction of Pq gradually 
increases as the vortex-vortex interactions dominate at 
higher fields. At the matching fields, the random array 
of pins has no peaks or enhancements in M(H) or Pk- 
This suggests that the presence of peaks in the periodic 
pinning arrays are due to the commensurability effects 
with the pinning substrate. 

The maximum value of the absolute value of M{H), 
\M(H)\, is 0.0095 $o/A 2 for the random pinning in Fig. 3, 
while it is 0.015 $0/ A 2 for the triangular array and square 
array. The latter value is about 1.5 times larger than that 
found for the random pinning case. This enhancement 
of M(H) occurs only for a limited range of fields. The 
M{H) for the triangular pinning array falls to the same 
value as M(H) for the random pinning array at H sa 
0.30 $o/A 2 , which is less than 25^. The M(H) for the 
square pinning array remains higher than the M(H) of 
the random pinning array until H w 0A5 ^>q/ X 2 . This 
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higher value of H at which the drop occurs for the square 
pinning array is due to the strong commensurability at 
the 2nd MF for the square pinning array, whereas for the 
triangular array the second MF is a less stable defective 
honeycomb lattice. For the triangular pinning array at 
the third MF, M(3Bj,) » 0.004 $ /A 2 , while the random 
pinning gives M(3B <I> ) « 0.0025 $ /A 2 . For fields higher 
than the fourth MF, M(H) is of the same order for the 
three pinning array geometries studied here. 

VII. REMARKS ON FINITE-SIZE EFFECTS 

Regarding finite size effects, we would like to empha- 
size that we have conducted simulations in samples that 
vary in size from 36A x 36A up to 72A x 72A, and we ob- 
serve the same features in all our simulations regardless 
of the system size. We have also done simulations with 
different pinning strengths and observe the same peaks 
in M(H) and Pq(H). This reproducibility in the peaks 
in different simulations suggest that the peaks are not 
merely fluctuations but are robust and reproducible re- 
sults. To further address this issue we have included in 
Fig. 3 both M(H) and Pk(H) for a system with the same 
pinning parameters as in the first two plots of M(H) but 
with pinning placed randomly. In this plot no peaks are 
visible in M(H) beyond the initial peak nor are any peaks 
visible in Pk (H) . The same behavior for the random ar- 
ray is observed for different sized systems. If the peaks in 
M(H) in systems with square and triangular pinning are 
due to finite size effects such as commensurability with 
the boundary conditions, then peaks in M(H) and Pk for 
system with the same size and boundary conditions but 
with random pinning should appear as well. 

The absence of any peaks in M{H) and Pk for the sys- 
tem with random pinning strongly suggests that peaks in 
these quantities for the square and triangular pinning ar- 
ray are due to commensurability effects with the pinning 
lattice only. It is important to stress that in our simula- 
tions, our analytical results and experimentally observed 
vortex lattice (VL) configurations are all consistent with 
each other. 

VIII. CONCLUSION 

To summarize, we have studied VLs interacting with 
periodic pinning arrays in which interstitial pinning is rel- 
evant above the first MF. We have shown that this system 
behaves considerably differently from atoms on surfaces 
or Josephson-junction arrays. A rich variety of distinct 
VLs can be stabilized including several novel partially- 
ordered lattices. We have also derived commensurability 
conditions for MFs at which stable ordered VLs appear. 
For the triangular pinning array these commensurability 



conditions are in excellent agreement with our simula- 
tions, while for the square array the commensurability 
conditions work for low fields (up to the 10th match- 
ing field). Our simulations are in excellent agreement 
with recent imaging experiments [Q and are robust over 
a wide range of parameters and system sizes. Our predic- 
tions can be tested with Lorentz microscopy techniques, 
Bitter decoration techniques, and by looking for a spe- 
cific sequence of peaks in magnetization measurements. 
These phases should be accessible to other systems with 
periodic pinning, including charged colloidal particles in 
a periodic array of optical traps || , and magnetic bubble 
arrays interacting with patterned substrates. 
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